version 7
#delimit;
set more off;
log using boehmke2006pa-analyze.log, text replace;

/*	************************************************************************	*/
/*     	File Name:	boehmke2006pa-analyze.do					*/
/*     	Date:   	May 12, 2006							*/
/*      Author: 	Frederick J. Boehmke						*/
/*      Purpose:	This file runs the analysis of the position timing and content 	*/
/*			data. There are separate models as well as estimation of the	*/
/*			four SUDCD models.						*/
/*      Input File:	boehmke2006pa.dta						*/
/*      Output File:	boehmke2006pa-analyze.log 					*/
/*	Version:	Stata 7 or above.				 		*/
/*	************************************************************************	*/


	/********************************************************/
	/* Below are the likelihood functions for the various 	*/
	/* estimators used here. 				*/
	/* 							*/
	/* 1. expdisc: This is the discrete exponential model.	*/
	/* 2. expwbl:  This is the Weibull SUDCD estimator.	*/
	/* 3. prolgn:  This is the log-normal SUDCD estimator.	*/
	/* 							*/
	/* Note that for the latter two models, you must 	*/
	/* specify the _rtcens variable	(as on line 140 below), */
	/* since it is hard-coded into the likelihoods. 	*/
	/********************************************************/


program define expdisc;

	version 7;
	args lnf theta1;
	quietly replace `lnf' = ln(exp(-exp(-`theta1')))	if $ML_y1==1;
	quietly replace `lnf' = ln(1-exp(-exp(-`theta1')))	if $ML_y1==0;

end;


program define expwbl;

	version 7;
	args lnf theta1 theta2 theta3 theta4;

	tempvar lambda1 lambda2 alphastar durdep P0_cond P_marg P_cum P_cens0 P_cens1;

	quietly gen double `lambda1'	= exp(-`theta1');
	quietly gen double `lambda2'	= exp(-`theta2');
	quietly gen double `alphastar' 	= (exp(2*`theta3')-1)/(exp(2*`theta3')+1);
	quietly gen double `durdep' 	= exp(`theta4');

	quietly gen double `P0_cond'	= 1 - exp(-`lambda1')*(1+`alphastar'*(2*exp(-($ML_y2*`lambda2')^`durdep')-1)*(exp(-`lambda1')-1));
	quietly gen double `P_marg'	= `lambda2'*`durdep'*((`lambda2'*$ML_y2)^(`durdep'-1))*exp(-(`lambda2'*$ML_y2)^`durdep');
	quietly gen double `P_cum'	= (1-exp(-`lambda1'))*(1-exp(-($ML_y2*`lambda2')^`durdep'))*(1+`alphastar'*exp(-($ML_y2*`lambda2')^`durdep'-`lambda1'));
	quietly gen double `P_cens0'	= (1-exp(-`lambda1')) - `P_cum';
	quietly gen double `P_cens1'	= 1 - (1-exp(-`lambda1')) - (1-exp(-($ML_y2*`lambda2')^`durdep')) + `P_cum';

	quietly replace `lnf' = ln(`P0_cond')   + ln(`P_marg') 	if $ML_y1==0 & _rtcens==0;
	quietly replace `lnf' = ln(1-`P0_cond') + ln(`P_marg') 	if $ML_y1==1 & _rtcens==0;
	quietly replace `lnf' = ln(`P_cens0') 			if $ML_y1==0 & _rtcens==1;
	quietly replace `lnf' = ln(`P_cens1') 			if $ML_y1==1 & _rtcens==1;

end;


program define prolgn;

	version 7;
	args lnf theta1 theta2 theta3 theta4;

	tempvar lambda2 rhostar sigma P0_cond P_marg P_cens0 P_cens1;

	quietly gen double `lambda2'	= exp(-`theta2');
	quietly gen double `rhostar' 	= (exp(2*`theta3')-1)/(exp(2*`theta3')+1);
	quietly gen double `sigma' 	= exp(`theta4');

	quietly gen double `P0_cond'	= 1 - norm((`theta1'+ln($ML_y2*`lambda2')*`rhostar'/`sigma')/sqrt(1-`rhostar'^2));
	quietly gen double `P_marg'	= normd(ln($ML_y2*`lambda2')/`sigma')/`sigma';
	quietly gen double `P_cens0'	= binorm( -ln($ML_y2*`lambda2')/`sigma' , -`theta1', -`rhostar' );
	quietly gen double `P_cens1'	= binorm( -ln($ML_y2*`lambda2')/`sigma' ,  `theta1',  `rhostar' );

	quietly replace `lnf' = ln(`P0_cond')   + ln(`P_marg') 	if $ML_y1==0 & _rtcens==0;
	quietly replace `lnf' = ln(1-`P0_cond') + ln(`P_marg') 	if $ML_y1==1 & _rtcens==0;
	quietly replace `lnf' = ln(`P_cens0') 			if $ML_y1==0 & _rtcens==1;
	quietly replace `lnf' = ln(`P_cens1') 			if $ML_y1==1 & _rtcens==1;

end;


	/********************************************************/
	/* Run the separate vote and duration models first. 	*/
	/* These models treat failures on day of vote as right- */
	/* censored (using the position indicator). Also 	*/
	/* tabulate predicted and observed values for PRE. 	*/
	/********************************************************/


use boehmke2006pa;

  probit vote contdiff mexbordr pscenter hhcenter partyid numdiff;

  	predict yhat_prob if e(sample), p;
  	recode 	yhat_prob 0/0.5=0 0.5/1=1;
  	tab 	yhat_prob vote, matcell(crosstab);


  ml model lf expdisc (vote = contdiff mexbordr pscenter hhcenter partyid numdiff);

	ml search;
	ml maximize;

	predict ystar_exp if e(sample), xb;
	generat yhat_exp = exp(-exp(-ystar_exp));
  	recode 	yhat_exp 0/0.5=0 0.5/1=1;

	tab yhat_exp vote, matcell(crosstab);


  stset timing, failure(position);

  streg corptpct labtpct mexbordr dleader rleader ncomact ideol pscenter inter1 hhcenter inter2, 
	dist(weibull) time;

  streg corptpct labtpct mexbordr dleader rleader ncomact ideol pscenter inter1 hhcenter inter2, 
	dist(lognormal) time;



	/********************************************************/
	/* Run the four SUDCD models. Note that right-censoring */
	/* is hard coded into the likelihood functions, which 	*/
	/* requires explicitly declaring the _rtcens variable.	*/
	/* To run without right-censoring, just set this 	*/
	/* variable equal to zero.				*/
	/********************************************************/

generat _rtcens = rtcensr;

	/********************************************************/
	/* First the two with no correlation parameterization. 	*/
	/********************************************************/


ml model lf expwbl (vote= contdiff mexbordr pscenter hhcenter partyid numdiff) 
  (timing=corptpct labtpct mexbordr dleader rleader ncomact ideol pscenter inter1 hhcenter inter2) 
  /alphastar /durdep;

		/********************************************************/
		/* Initial parameter search may send alphastar to 	*/
		/* boundary, so push it back into reasonable values so	*/
		/* that the maximization won't crash.			*/
		/********************************************************/

	ml search;
	  ml init /alphastar=0;
	ml maximize, difficult;

		/********************************************************/
		/* The _diparm commands convert the estimation results 	*/
		/* for the correlation parameter into its bounded	*/
		/* range (-1<alpha<1 for weibull SUDCD with rho=alpha/4 */
		/* and -1<rho<1 for log-normal SUDCD).		 	*/
		/********************************************************/

	_diparm alphastar, label("alpha") function((exp(2*(@))-1)/(exp(2*(@))+1))      
	  derivative(4*exp(2*@)/(exp(2*@)+1)^2);

	_diparm alphastar, label("rho")   function(((exp(2*(@))-1)/(exp(2*(@))+1))/4)  
	  derivative((4*exp(2*@)/(exp(2*@)+1)^2)/4);


ml model lf prolgn (vote= contdiff mexbordr pscenter hhcenter partyid numdiff) 
  (timing=corptpct labtpct mexbordr dleader rleader ncomact ideol pscenter inter1 hhcenter inter2) 
  /rhostar /sigma;

	ml search;
	ml maximize, difficult;

	_diparm rhostar, label("rho") function((exp(2*(@))-1)/(exp(2*(@))+1))  
	  derivative(4*exp(2*@)/(exp(2*@)+1)^2);


	/********************************************************/
	/* Now the two with the correlation parameterization. 	*/
	/********************************************************/


ml model lf expwbl (vote= contdiff mexbordr pscenter hhcenter partyid numdiff) 
  (timing=corptpct labtpct mexbordr dleader rleader ncomact ideol pscenter inter1 hhcenter inter2) 
  (alphastar: partyid) /durdep;

		/********************************************************/
		/* As before, put constant and partyid parameters back 	*/
		/* to reasonable values after initial parameter search 	*/
		/* so that the maximization won't crash.		*/
		/********************************************************/

	ml search;
	  ml init alphastar:_cons=0;
	  ml init alphastar:partyid=0;
	ml maximize, difficult;

		/********************************************************/
		/* Test for significant correlation among Democrats.	*/
		/********************************************************/


	display "Correlation for Republicans (rho): "
		((exp(2*([alphastar]_b[_cons]))-1)
		/(exp(2*([alphastar]_b[_cons]))+1))/4;

	display "Correlation for Democrats (rho): "   	
		((exp(2*([alphastar]_b[_cons] + [alphastar]_b[partyid]))-1)
		/(exp(2*([alphastar]_b[_cons] + [alphastar]_b[partyid]))+1))/4;

	test [alphastar]partyid + [alphastar]_cons = 0;


ml model lf prolgn (vote= contdiff mexbordr pscenter hhcenter partyid numdiff) 
  (timing=corptpct labtpct mexbordr dleader rleader ncomact ideol pscenter inter1 hhcenter inter2) 
  (rhostar: partyid) /sigma;

	ml search;
	ml maximize, difficult;

	_diparm rhostar, label("rho") function((exp(2*(@))-1)/(exp(2*(@))+1)) 
	  derivative(4*exp(2*@)/(exp(2*@)+1)^2);

	display "Correlation for Republicans (rho): "
		((exp(2*([rhostar]_b[_cons]))-1)
		/(exp(2*([rhostar]_b[_cons]))+1));

	display "Correlation for Democrats (rho): "   	
		((exp(2*([rhostar]_b[_cons] + [rhostar]_b[partyid]))-1)
		/(exp(2*([rhostar]_b[_cons] + [rhostar]_b[partyid]))+1));

	test [rhostar]partyid + [rhostar]_cons = 0;


log close;
clear;
exit;